rifle = function(data, unit, time, leader, outcome, nsim){

    set.seed(21024)
  
    rifle_run = function(data, unit, time, leader, outcome){
      
            options(dplyr.summarise.inform = FALSE)
    
            outcome_data = as.data.frame(data) %>%  
            select(unit, time, outcome) 
            colnames(outcome_data)[1:3] = c('unit', 'time', 'outcome')
                
            outcome_data = outcome_data %>%  
              group_by(unit, time) %>%  
              summarise(outcome = mean(outcome))
                
            outcome_data = outcome_data %>%  
              mutate(time = time - (min(time))+1 )
                
            shuffle_leader = function(data, unit, time, leader){
              
              data = data %>%
                select(unit, time, leader)
              colnames(data)[1:3] = c('unit', 'time', 'leader')
              
              start_year_true = min(data$time)
              
              data_shuffle = data %>% 
                group_by(unit, leader) %>%  
                mutate(min_year = min(time)) %>%  
                group_by(unit) %>%  
                mutate(leader_order = rank(min_year, ties.method = 'min')) %>%  
                mutate(order = dense_rank(leader_order)) %>%  
                group_by(leader, unit) %>%
                mutate(tenure = n()) %>%  
                group_by(leader, unit) %>%  
                summarise(tenure = mean(tenure), 
                          order = mean(order)) %>%
                group_by(unit) %>%  
                mutate(
                  shuffle_order = sample(n(), n(), replace = F)
                ) %>%  
                group_by(unit) %>%  
                mutate(lag_tenure = lag(tenure, 1, order_by = shuffle_order))  %>%  
                mutate(start_time = if_else(shuffle_order==1, 1, 
                                            1+ lag_tenure)) 
              
              data_shuffle = data_shuffle %>%  
                arrange(unit, shuffle_order) %>%  
                group_by(unit) %>% 
                mutate(first_year = ifelse(shuffle_order == 1, 1, cumsum(replace_na(lag_tenure, 0))+1), 
                       last_year = (first_year + tenure)-1)
              
              # Create a sequence of years between the first year and the last year
              year_sequence <- map2(
                data_shuffle$first_year,
                data_shuffle$last_year,
                ~ seq(.x, .y)
              )
              
              # Expand the original dataframe to include all the years each president served
              expanded_df <- tibble(
                name = rep(data_shuffle$leader, lengths(year_sequence)),
                unit = rep(data_shuffle$unit, lengths(year_sequence)),
                time = unlist(year_sequence)
              )
              
              
              return(expanded_df)
              
            }
            
            
    shuffle_data = shuffle_leader(data, unit, time, leader)
    
    shuffle_data_run = merge(shuffle_data, outcome_data,
                             by = c('unit', 'time'), all.y = T)
    
    return(summary(lm(outcome ~ factor(name), data = shuffle_data_run))$r.sq)
    
    }

            map_progress <- function(.x, .f, ..., .id = NULL) {
              .f <- purrr::as_mapper(.f, ...)
              pb <- progress::progress_bar$new(total = length(.x), format = " [:bar] :current/:total (:percent) eta: :eta", force = TRUE)
              
              f <- function(...) {
                pb$tick()
                .f(...)
              }
              purrr::map(.x, f, ..., .id = .id)
            }
result_list = unlist(map_progress(seq_len(nsim), ~rifle_run(data, unit, time, leader, outcome), 
              progress = TRUE))


outcome_data = as.data.frame(data) %>%  
  select(unit, time, outcome, leader) 
colnames(outcome_data)[1:4] = c('unit', 'time', 'outcome', 'leader')
outcome_data = outcome_data %>%  
  group_by(unit, time, leader) %>%  
  summarise(outcome = mean(outcome))
outcome_data = outcome_data %>%  
  mutate(time = time - (min(time))+1 )

real_r2 = summary(lm(outcome ~ factor(leader), data = outcome_data))$r.sq

p_value = (sum(as.numeric(result_list>real_r2)))/nsim

return(c(result_list, p_value))

}
